For burn projects that were wrapped up in less than a week:
a.) How well did the number of days for which thermal anomalies were detected corroborate with the number of days that burners indicated in their project reporting?
b.) Is there an apparent detection limit (e.g., acres burned or some other metric) under which the satellite does not seem to perform well in terms of picking up thermal anomaly? How does the finding different for broadcast burn vs. pile burn projects? How does the apparent detection limit compare to the US EPA’s Rx burn reporting thresholds (50-acre for broadcast burn; 25-acre for pile burn)?
# VIIRS fire detects
viirs <- read_rds(file.path(ref_path, 'VIIRS/CA_2020/viirs_nonAg.rds'))
# CalMAPPER data
CM <- read_rds(file.path(ref_path, 'CalMAPPER/treatment-activities.rds'))
Just get burn data out of CalMAPPER.
# check if cultural burning can be either pile burn or broadcast. all broadcast this year at least.
CM %>%
as_tibble %>%
filter(ACTIVITY_DESCRIPTION == 'Cultural Burning') %>%
distinct(ACTIVITY_DESCRIPTION, ACTIVITY_NAME)
## # A tibble: 2 × 2
## ACTIVITY_DESCRIPTION ACTIVITY_NAME
## <chr> <chr>
## 1 Cultural Burning Prescribed Fire - Cultural Burning
## 2 Cultural Burning Prescribed Fire - Oak Woodland Restoration
# CalMAPPER -- <= 7 days, trts include burns
CM_burn <- CM %>%
# create columns for burn activity
mutate(burn = ACTIVITY_DESCRIPTION %in%
c('Broadcast Burn', 'Cultural Burning', 'Pile Burning'),
pile_burn = ACTIVITY_DESCRIPTION == 'Pile Burning',
broadcast_burn = ACTIVITY_DESCRIPTION %in%
c('Broadcast Burn', 'Cultural Burning'),
.after = AQ_ID) %>%
# get duration of burn activities. +1 to duration (e.g. 2/1-2/1 = 0. It should be 1 day.)
mutate(duration = difftime(ACTIVITY_END, ACTIVITY_START, units = 'days') + 1) %>%
# include only burn activity
filter(burn == T) %>%
# get area of treatment too. might be useful
mutate(trt_acres = units::set_units(st_area(.), 'acres') ) %>%
# keep only the useful columns
select(TREATMENT_ID, AQ_ID, burn:broadcast_burn, TREATMENT_NAME, ACTIVITY_NAME, ACTIVITY_DESCRIPTION, ACTIVITY_START, ACTIVITY_END, duration, ACTIVITY_STATUS, TREATED_ACRES, trt_acres, BROAD_VEG_TYPE, AGENCY_NAME, FY)
Just look at projects from 2020. CalMAPPER is supposed to be complete only FY19/20 and up. I downloaded VIIRS data from Jan1-Dec31. Check out the projects that span calendar years in and out of 2020… If project spans calendar year, the start date is always on Dec 31st. Seems like a default value.
Update: Dec 31 1600 PST is the same as Jan 1 0000 UTC. It seems like a default value whenever the start date isn’t known. So, remove the 2021 start dates, but 12/31/2019 starts are ok.
# projects spanning 2020 calendar year
CM_burn %>%
filter(year(ACTIVITY_END) > 2020 & year(ACTIVITY_START) == 2020 |
year(ACTIVITY_START) < 2020 & year(ACTIVITY_END) == 2020) %>%
pull(ACTIVITY_START) %>%
unique
## [1] "2019-12-31 16:00:00 PST" "2020-12-31 16:00:00 PST"
# filter just 2020
CM_burn2020 <- CM_burn %>%
filter(year(ACTIVITY_END) == 2020)
Let’s also just look at the distribution of start and end dates. Are their common dates that the records tend to show? Spikes in the figure would indicate inaccurate dates.
CalMapper_dates <- cowplot::plot_grid(
ggplot(CM_burn, aes(ACTIVITY_START)) +
geom_bar(),
ggplot(CM_burn, aes(ACTIVITY_END)) +
geom_bar(),
nrow = 2
)
CalMapper_dates
ggsave(plot = CalMapper_dates,
filename = 'tmp_outputs/CalMAPPER_activity_dates.png', width = 6, height = 4)
This data set has 2 treatments that have the exact same geometries, even though they have different ID numbers. This will cause issues when trying to count the number of unique fire detect days within the polygons. e.g. one treatment has a burn in january and the other treatment has burns from feb-may. The additional detects will seem anomolous for the first treatment, when in fact it should be part of the second treatment.
Since there are only 2 treatments that are duplicated, manually merge them. In the future, manual cleaning up of the data should be avoided.
# calmapper treatment ids.
CM_trtIDs <- CM_burn2020 %>%
group_by(TREATMENT_ID) %>%
summarise()
# find which rows have duplicates
mat_eq <- st_equals(CM_trtIDs, sparse = F, remove_self = T)
duplicated_trts <- CM_trtIDs$TREATMENT_ID[apply(mat_eq, 1, sum) == 1]
CM_burn2020 %>% filter(TREATMENT_ID %in% c(duplicated_trts))
## Simple feature collection with 22 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13463530 ymin: 4732515 xmax: -13461520 ymax: 4737232
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
## TREATMENT_ID AQ_ID burn pile_burn broadcast_burn TREATMENT_NAME
## 1 9553 60497 TRUE TRUE FALSE N-05-19 Phase I
## 2 14544 62420 TRUE TRUE FALSE North Fork American River
## 3 14544 63060 TRUE TRUE FALSE North Fork American River
## 4 14544 67008 TRUE TRUE FALSE North Fork American River
## 5 14544 68189 TRUE TRUE FALSE North Fork American River
## 6 14544 74051 TRUE TRUE FALSE North Fork American River
## 7 14544 74066 TRUE TRUE FALSE North Fork American River
## 8 14544 74898 TRUE TRUE FALSE North Fork American River
## 9 14544 77458 TRUE TRUE FALSE North Fork American River
## 10 14544 78102 TRUE TRUE FALSE North Fork American River
## ACTIVITY_NAME ACTIVITY_DESCRIPTION ACTIVITY_START
## 1 2019 CNG & Fuels Crew Pile Burning Pile Burning 2020-01-05 16:00:00
## 2 2020 CNG Pile Burning Pile Burning 2019-12-31 16:00:00
## 3 2020 CNG Pile Burning Pile Burning 2020-01-19 16:00:00
## 4 2020 CNG Pile Burning Pile Burning 2020-01-26 16:00:00
## 5 2020 CNG Pile Burning Pile Burning 2020-02-02 16:00:00
## 6 2020 CNG Pile Burning Pile Burning 2020-02-09 16:00:00
## 7 2020 CNG Pile Burning Pile Burning 2020-02-23 16:00:00
## 8 2020 CNG Pile Burning Pile Burning 2020-03-01 16:00:00
## 9 2020 CNG Pile Burning Pile Burning 2020-03-08 17:00:00
## 10 2020 CNG Pile Burning Pile Burning 2020-03-15 17:00:00
## ACTIVITY_END duration ACTIVITY_STATUS TREATED_ACRES trt_acres
## 1 2020-01-09 16:00:00 5 days Complete 1.5 1501.376 [acres]
## 2 2020-01-16 16:00:00 17 days Complete 6.5 1501.376 [acres]
## 3 2020-01-23 16:00:00 5 days Complete 7.0 1501.376 [acres]
## 4 2020-01-30 16:00:00 5 days Complete 6.5 1501.376 [acres]
## 5 2020-02-06 16:00:00 5 days Complete 6.0 1501.376 [acres]
## 6 2020-02-20 16:00:00 12 days Complete 8.0 1501.376 [acres]
## 7 2020-02-27 16:00:00 5 days Complete 3.0 1501.376 [acres]
## 8 2020-03-05 16:00:00 5 days Complete 4.0 1501.376 [acres]
## 9 2020-03-12 17:00:00 5 days Complete 4.0 1501.376 [acres]
## 10 2020-03-19 17:00:00 5 days Complete 1.0 1501.376 [acres]
## BROAD_VEG_TYPE AGENCY_NAME FY SHAPE
## 1 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 2 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 3 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 4 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 5 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 6 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 7 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 8 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 9 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 10 Timber CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
# change the duplicated treatment ID
CM_burn2020$TREATMENT_ID[CM_burn2020$TREATMENT_ID == 9553] <- 14544
# rerun to update
CM_trtIDs <- CM_burn2020 %>%
group_by(TREATMENT_ID) %>%
summarise()
Export this data so that you can use it for the other questions.
write_rds(CM_burn2020, file.path(ref_path, 'CalMAPPER/burnactivities_2020.rds'))
write_rds(CM_trtIDs, file.path(ref_path, 'CalMAPPER/treatmentgeometries_2020.rds'))
Since I’m going to filter the VIIRS detections based on the CalMAPPER polygons, I should create a dataset with polygons with an extra buffer, where the radius is half that of the VIIRS resolution. This is just in case there’s a VIIRS detection right outside the polygon–I don’t want to miss that. And export.
# buffer polygons
st_crs(CM_burn2020)$units # make sure buffer will be in meters
## [1] "m"
CM_burn2020_buffered <- st_buffer(CM_burn2020, 375/2)
CM_trtIDs_buffered <- CM_burn2020_buffered %>%
group_by(TREATMENT_ID) %>%
summarise()
# export
write_rds(CM_burn2020_buffered, file.path(ref_path, 'CalMAPPER/burnactivities_2020_buffered.rds'))
write_rds(CM_trtIDs_buffered, file.path(ref_path, 'CalMAPPER/treatmentgeometries_2020_buffered.rds'))
Finally, keep treatments where burn activity was <= 7 days.
CM_7days <- CM_burn2020_buffered %>%
filter(duration <= 7)
Filter out VIIRS data outside of the the CalMAPPER polygons. I’m using the buffered CalMAPPER polygons.
# get the viirs points contained in the treatment polygons
viirs_t <- st_transform(viirs, st_crs(CM_7days))
# calmapper treatment ids.
CM_7days_trtIDs <- filter(CM_trtIDs_buffered, TREATMENT_ID %in% unique(CM_7days$TREATMENT_ID))
# CM_7days_trtIDs_unbuff <- filter(CM_trtIDs, TREATMENT_ID %in% unique(CM_7days$TREATMENT_ID))
# mapview(CM_7days_trtIDs) + mapview(CM_7days_trtIDs_unbuff)
# find those that intersect
viirs_int <- viirs_t %>%
filter(CONFIDENCE != 'l') %>%
st_intersection(CM_burn2020_buffered) %>%
mutate(rowname = rownames(.),
pointID = floor(as.numeric(rowname)),
overlap_trts = grepl('\\.', rowname),
.before = LATITUDE) %>%
group_by(pointID) %>%
mutate(overlap_trts = any(overlap_trts == T)) %>%
ungroup()
Find if any points sit on overlapping treatment. It looks like a handful of VIIRS points intersect overlapping polygons. This seems to often happen when there’s burn prep ahead of a broadcast burn. Out of the 1463 VIIRS points that intersect the CalMAPPER polygons, only 101 (7% of detects) are over overlapping treatments. Ignore for now and add a caveat for potential artifacts?
# find if any points sit on overlapping treatment
viirs_overlap <- filter(viirs_int, overlap_trts == T)
# number of viirs points vs those that overlap polygons
viirs_int %>%
as_tibble() %>%
distinct(pointID, .keep_all = T) %>%
count(overlap_trts)
## # A tibble: 2 × 2
## overlap_trts n
## <lgl> <int>
## 1 FALSE 721
## 2 TRUE 1404
# visualize
mapview::mapview(viirs_overlap, cex = 3, zcol = 'TREATMENT_ID', legend = F) +
mapview::mapview(CM_7days, zcol = 'TREATMENT_NAME', legend= F, alpha.regions = 0.3)